import pandas as pd
import numpy as np
from pandas import read_csv
from matplotlib import pyplot
import matplotlib.pyplot as plt
Задание 9:¶
Для заданной страны необходимо построить прогнозы на ближайшие 5 лет показателей деловой активности (индекс промышленного производства, индекс деловой активности, ВВП в постоянных ценах, уровень безработицы). Использовать месячные данные. В качестве возможных факторов использовать: фондовый индекс, уровень процентных ставок, отраслевые индексы цен, курс национальной валюты к доллару США, уровень цен на основные сырьевые товары (нефть и золото), дефлятор ВВП, индекс потребительских цен и основные денежные агрегаты (денежная база, денежная масса (M1, M2)).
Исследовать, улучшают ли выбранные факторы прогноз показателей деловой активности в сравнении с одномерной моделью временного ряда?
Для нашего исследования деловой активности мы решили использовать следующие данные:
Какие данные используются?¶
- Индекс деловой активности в производственном секторе (PMI) США от ISM
- Композитный индекс деловой активности (PMI) в США от Markit
- Фьючерс на нефть Brent
- Фьючерс на золото
- Уровень безработицы
- U.S. Gross Domestic Product
- S&P 500 (SPX)
- Объем производства %
- Курс Доллара к Евро
Среди показателей деловой активности мы выбрали два похожих по содержанию показаткелей - индексы деловой активности в производственном секторе и композитный индекс деловой активности.
Расскажем немного о них:¶
PMI — это макроэкономический показатель, позволяющий судить о стабильности экономики. Во всем мире для оценки состояния экономики используют индексы деловой активности. Наибольшее распространение они получили в США, Китае и странах Еврозоны. Эти показатели имеют большое влияние на экономические процессы, но главным образом на рынок ценных бумаг.
Индекс PMI рассчитывается на основании опроса топ-менеджеров ключевых производственных предприятий. Ответы даются на основании сравнения существующих показателей со статистикой истекшего периода (календарный месяц). Значение PMI, которое представляет собой число от 0 до 100, определяет собой оценку изменения Валютного курса страны, Темпы роста ВВП, стоимость ценных бумаг и других ведущих показателей.
Разновидности индекса PMI¶
- Активность в производственном секторе
- Активность в сфере услуг
- Активность в непроизводственном секторе
- Композитный индекс
В своем проекте мы решили взять для анализа и предсказания два наиболее популярных индекса- композитный и произвеодственный в качестве предсказываемых переменных.Именно эти индексы наиболее чуствительны к экономическим колебаниям и изменениям в экономике стран. Мы также выбрали в качестве объясныющих переменных основыной индекс S&P 500, который включает в себя акции наиболее крупных Американских компаний, большинство менеджеров которых и определяют индекс PMI. Мы также собрали данные по ВВП страны и основные котировки на сырьевом рынке и рынке драг.металлов: Brent и Золото, а также ИПЦ США.
Цель исследования :¶
В кругах финансовых аналитиков ходят разные мнения на счет направления влияния индекса деловой активности и определяющих факторов этого индекса.Кто-то говорит, что оценка менеджерами "перспектив" влияет на ВВП, акции и иные финансовые индексы, кто-то считает наоборот.В данном проекте мы попытаемся проанализировать временные ряды, определить их свойства и построить прогноз на ближайшие 5 лет индексов PMI.
Загрузка и визуализация данных:¶
import matplotlib.pyplot as plt
import plotly.express as px
from datetime import datetime
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
df['month_year'] = pd.to_datetime(df['Date']).dt.to_period('M')
#df=df.sort_values(by='Date')
fig = px.line(df, x='Date', y='PMIcomp')
fig.add_hline(y=50.0,line_dash="dash", line_color="red")
fig.show()
fig = px.line(df, x='Date', y='Vol')
fig.show()
fig = px.line(df, x='Date', y='PMItot')
fig.add_hline(y=50.0,line_dash="dash", line_color="red")
fig.show()
fig = px.line(df, x='Date', y='SP')
fig.show()
fig = px.line(df, x='Date', y='GDP')
fig.show()
fig = px.line(df, x='Date', y='Gold')
fig.show()
fig = px.line(df, x='Date', y='Brent')
fig.show()
fig = px.line(df, x='Date', y='Unemployment')
fig.show()
fig = px.line(df, x='Date', y='USD-EUR')
fig.show()
#df.PMIcomp = df.PMIcomp.astype(int)
#df.PMItot= df.PMItot.astype(int)
# Everything in one
fig = px.line(df, x='Date', y=['Vol','USD-EUR','SP','GDP','Brent','Gold','PMItot','PMIcomp'])
fig.show()
Анализ на стационарность
from statsmodels.tsa.arima.model import ARIMA
import statsmodels.tsa.stattools
import statsmodels.api as sm
from statsmodels.tsa.stattools import adfuller
alpha=0.05
#perform augmented Dickey-Fuller test
result = adfuller(df['SP'], autolag='AIC')
print(f'ADF Statistic: {result[0]}')
print(f'p-value: {result[1]}')
for key, value in result[4].items():
print('Critial Values:')
print(f' {key}, {value}')
if result[1] < 0.05:
print('Ряд стационарен')
else:
print('Ряд не стационарен')
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['GDP'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Gold'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Brent'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['PMItot'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['PMIcomp'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['USD-EUR'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
KPSS -test. Так как индексы деловой активности строго ограницены [0-100] и колеблются около 50 - мы должны провести иной тест на стационарность, диагностирующий наличие либо отсутствие тренда.
from statsmodels.tsa.stattools import kpss
import warnings
warnings.filterwarnings("ignore")
def kpss_test(series, **kw):
statistic, p_value, n_lags, critical_values = kpss(series, **kw)
# Format Output
print(f'KPSS Statistic: {statistic}')
print(f'p-value: {p_value}')
print(f'num lags: {n_lags}')
print('Critial Values:')
for key, value in critical_values.items():
print(f' {key} : {value}')
print(f'Result: The series is {"not " if p_value < 0.05 else ""}stationary')
kpss_test(df['PMItot'])
kpss_test(df['PMIcomp'])
kpss_test(df['Brent'])
kpss_test(df['Gold'])
kpss_test(df['SP'])
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['Unemployment']]
ts1['Unemployment'] =ts1['Unemployment'].str.replace('%', '')
ts1=pd.to_numeric(ts1['Unemployment'])
kpss_test(ts1)
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(ts1)
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
Выводы:¶
Мы проанализировали стационарность и наличие тренда в изучаемых временных рядах. Наши индексы деловой активности явялются стационарными рядами с детерминированным трендом, GDP по тесту ADF также является стационарным детерминированным временным рядом. Особое внимание стоит обратить на уровень безработицы, который имеет нестационарный детерминированный тренд. Индекс S&P 500, фьючерсы на нефть и золото явялются нестационарными рядами со стохастическим временным трендом.
На графике ВВП мы видим характерные черты сезонности, проверим нашу гипотезу. Обнаружение и удаление сезонности.
Есть два способа:
1 Дифференциация: разница в значении во время определенного количества лагов.
2 Декомпозиция: моделирование трендов и сезонности отдельно удаляет их
y(t) = Level + Trend + Seasonality + Noise
from statsmodels.tsa.seasonal import seasonal_decompose
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
df['Date'] = pd.to_datetime(df['Date'])
df.set_index('Date', inplace=True)
df=df[['GDP']]
analysis = df[['GDP']].copy()
decompose_result_additive = seasonal_decompose(analysis, model="ad",period=12)
trend = decompose_result_additive.trend
seasonal = decompose_result_additive.seasonal
residual = decompose_result_additive.resid
decompose_result_additive .plot()
Как мы видим из графика сезонность есть.В дальнейшем мы будем использовать этот факт для составления прогноза для индексов деловой активности.Проверим наши исследования:
Подбор параметров ARMA- модели.Первые прогнозы:¶
ARIMA models can be applied only in stationary data.
from statsmodels.tsa.arima_model import ARIMA
import pmdarima as pm
from itertools import product
import pmdarima as pm
from statsmodels.graphics.tsaplots import acf, pacf, plot_acf, plot_pacf
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts=df[['GDP']]
ts['GDP'] = (ts['GDP'] - ts['GDP'].mean()) / ts['GDP'].std() #приведение к нормальному распределению
n_sample = len(ts['GDP'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts['GDP'].iloc[:pivot], ts['GDP'].iloc[pivot:]
arima_model = pm.auto_arima(x_train,start_p=1, start_q=1,
test='adf',
max_p=3, max_q=3, m=12,
start_P=0, seasonal=True,
d=None, D=1, trace=True,
error_action='ignore',
suppress_warnings=True,
stepwise=True)
arima_model.summary()
plot_pacf(ts).show()
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['Gold']]
ts1['Gold'] = (ts1['Gold'] - ts1['Gold'].mean()) / ts1['Gold'].std() #приведение к нормальному распределению
n_sample = len(ts1['Gold'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['Gold'].iloc[:pivot], ts1['Gold'].iloc[pivot:]
arima_model = pm.auto_arima(x_train,
start_p=1,d =0,start_q=0,
test="adf", supress_warnings = True,
trace=True)
arima_model.summary()
plot_pacf(ts1).show()
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['PMItot']]
ts1['PMItot'] = (ts1['PMItot'] - ts1['PMItot'].mean()) / ts1['PMItot'].std() #приведение к нормальному распределению
n_sample = len(ts1['PMItot'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['PMItot'].iloc[:pivot], ts1['PMItot'].iloc[pivot:]
arima_model = pm.auto_arima(x_train,
start_p=1,d =0,start_q=0,
test="adf", supress_warnings = True,
trace=True)
arima_model.summary()
plot_pacf(ts1).show()
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['PMIcomp']]
ts1['PMIcomp'] = (ts1['PMIcomp'] - ts1['PMIcomp'].mean()) / ts1['PMIcomp'].std() #приведение к нормальному распределению
n_sample = len(ts1['PMIcomp'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['PMIcomp'].iloc[:pivot], ts1['PMIcomp'].iloc[pivot:]
arima_model = pm.auto_arima(x_train,
start_p=1,d =0,start_q=0,
test="adf", supress_warnings = True,
trace=True)
arima_model.summary()
plot_pacf(ts1).show()
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['USD-EUR']]
ts1['USD-EUR'] = (ts1['USD-EUR'] - ts1['USD-EUR'].mean()) / ts1['USD-EUR'].std() #приведение к нормальному распределению
n_sample = len(ts1['USD-EUR'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['USD-EUR'].iloc[:pivot], ts1['USD-EUR'].iloc[pivot:]
arima_model = pm.auto_arima(x_train,
start_p=1,d =0,start_q=0,
test="adf", supress_warnings = True,
trace=True)
arima_model.summary()
plot_pacf(ts1).show()
#приведем данные к стационарному виду
df = pd.read_csv("final.csv", parse_dates=True,sep=';')
df['Unemployment'] = df['Unemployment'].str.replace('%', '')
df['Unemployment']=pd.to_numeric(df['Unemployment'])
df['Unemployment_dif']=(df['Unemployment']-df['Unemployment'].shift(1))
df[['Unemployment_dif']].dropna()
df= df[df['Unemployment_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Unemployment_dif'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
#приведем данные к стационарному виду
df['Brent_dif']=(df['Brent']-df['Brent'].shift(1)).dropna()
df[['Brent_dif']].dropna()
df= df[df['Brent_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Brent_dif'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
df['Gold_dif']=(df['Gold']-df['Gold'].shift(1)).dropna()
df[['Gold_dif']].dropna()
df= df[df['Gold_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['Gold_dif'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
df['SP_dif']=(df['SP']-df['SP'].shift(1)).dropna()
df[['SP_dif']].dropna()
df= df[df['SP_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['SP_dif'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
df['USD-EUR_dif']=(df['USD-EUR']-df['USD-EUR'].shift(1)).dropna()
df[['USD-EUR_dif']].dropna()
df= df[df['USD-EUR_dif'].notna()]
stat, pvalue, *_ = statsmodels.tsa.stattools.adfuller(df['USD-EUR_dif'])
stat, pvalue
if pvalue < alpha:
print(stat, pvalue,'Ряд стационарен')
else:
print(stat, pvalue,'Ряд не стационарен')
Первые прогнозы
import statsmodels.api as sm
import pmdarima as pm
from itertools import product
import pmdarima as pm
#df = pd.read_csv("final.csv", parse_dates=True,sep=';')
ts1=df[['PMItot']]
ts1['PMItot'] = (ts1['PMItot'] - ts1['PMItot'].mean()) / ts1['PMItot'].std() #приведение к нормальному распределению
n_sample = len(ts1['PMItot'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['PMItot'].iloc[:pivot], ts1['PMItot'].iloc[pivot:]
first_val = x_train.iloc[0]
x_train_diff = x_train.diff().dropna()
model = sm.tsa.arima.ARIMA(x_train_diff.values, order=(1,0,0)).fit()
x_pred_diff = model.forecast(len(x_test) - 1)
x_pred_diff = np.array([first_val] + x_pred_diff.tolist())
x_pred = pd.Series(np.cumsum(x_pred_diff))
plt.plot(np.arange(len(x_pred)) + len(x_train),x_test)
plt.plot(np.arange(len(x_pred)) + len(x_train), x_pred)
plt.plot(x_train)
plt.plot(x_train)
plt.legend(['y_test', 'y_test_pred', 'y_train','y_train_pred'])
ts1=df[['PMIcomp']]
ts1['PMIcomp'] = (ts1['PMIcomp'] - ts1['PMIcomp'].mean()) / ts1['PMIcomp'].std() #приведение к нормальному распределению
n_sample = len(ts1['PMIcomp'])
pivot = int(0.7 * n_sample)
x_train, x_test = ts1['PMIcomp'].iloc[:pivot], ts1['PMIcomp'].iloc[pivot:]
first_val = x_train.iloc[0]
x_train_diff = x_train.diff().dropna()
model = sm.tsa.arima.ARIMA(x_train_diff.values, order=(2,0,0)).fit()
x_pred_diff = model.forecast(len(x_test)-1 )
x_pred_diff = np.array([first_val] + x_pred_diff.tolist())
x_pred = pd.Series(np.cumsum(x_pred_diff))
plt.plot(np.arange(len(x_pred)) + len(x_train),x_test)
plt.plot(np.arange(len(x_pred)) + len(x_train), x_pred)
plt.plot(x_train)
plt.plot(x_train)
plt.legend(['y_test', 'y_test_pred', 'y_train','y_train_pred'])
Как мы видим, несмотря на то, что модели, подобранные алгоритмом, максимизируют точность - актуальный прогноз построить не получилось, поэтому расмотрим подробнее зависимые переменые:
Множественная регрессия временных рядов:
Общая регрессионная модель временных рядов с несколькими объясняющими переменными является расширением авторегрессионной модели с распределенными лагами на случай включения в число регрессоров запаздываний нескольких переменных. В общем случае регрессионная модель временных рядов расширяется на случай k дополнительных регрессов и включает q1:
- запаздывания первой объясняющей переменной, q1
- запаздывания второй объясняющей переменной и так далее.
Однако большинство наших временных рядов являются нестационнарными временными рядами.И прежде чем составлять какой-либо моножественный прогноз мы должны выяснить
1) Обладают ли ряды детерминированными или стохастическими трендами
2) Устранение нестационарности: "на практике для стационарности многомерного ряда достаточно стационарности составляющих его одномерных рядов."
3) Провести анализ на наличие структурных разрывов
4) Провести тест на причинность по Гренджеру
4) Проверить коинтегрированность индексов в случае отсутствия ложной регрессии между ними (тест ADF-GL/ тест Йенсена)
5) Построить VAR() модель : оценить параметры; Построить многомерный прогноз;
Cross Correlation Function (CCF) and ADF-GL test + Johansen test for co-integration
# CFF-test and Granger-Causality test
cff=sm.tsa.stattools.ccf(df['SP_dif'],df['PMItot'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['Gold_dif'],df['PMItot'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['Brent_dif'],df['PMItot'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['GDP'],df['PMItot'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['USD-EUR'],df['PMItot'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['Unemployment'],df['PMItot'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['SP_dif'],df['PMIcomp'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['Gold_dif'],df['PMIcomp'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['Brent_dif'],df['PMIcomp'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['GDP'],df['PMIcomp'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['Unemployment_dif'],df['PMIcomp'], adjusted=False)
cff[:7]
cff=sm.tsa.stattools.ccf(df['USD-EUR'],df['PMIcomp'], adjusted=False)
cff[:7]
Выводы:
Тест на причинность по Гренджеру T he p-value of each test is less than 5%, so it can be said that the any index is useful for forecasting Index of intepreners' activity .
from statsmodels.tsa.stattools import grangercausalitytests
grangercausalitytests(df[['PMItot', 'USD-EUR']], maxlag=[2])
grangercausalitytests(df[['PMItot', 'GDP']], maxlag=[1])
grangercausalitytests(df[['PMItot','Gold_dif']], maxlag=[1])
grangercausalitytests(df[['PMItot','Brent_dif']], maxlag=[1])
grangercausalitytests(df[['PMItot','SP_dif']], maxlag=[1])
grangercausalitytests(df[['PMIcomp','SP_dif']], maxlag=[1])
grangercausalitytests(df[['PMIcomp','Brent_dif']], maxlag=[1])
grangercausalitytests(df[['PMIcomp','GDP']], maxlag=[1])
grangercausalitytests(df[['PMIcomp','Unemployment_dif']], maxlag=[1])
grangercausalitytests(df[['PMIcomp','Gold_dif']], maxlag=[1])
grangercausalitytests(df[['PMIcomp','PMItot']], maxlag=[1])
Выводы
from statsmodels.tsa.vector_ar.vecm import coint_johansen
from statsmodels.tsa.vector_ar.vecm import coint_johansen
# dataframe of n series for cointegration analysis
import numpy as np
from statsmodels.tsa.tsatools import lagmat
def johansen(ts, lags):
"""
Calculate the Johansen Test for the given time series
"""
# Make sure we are working with arrays, convert if necessary
ts = np.asarray(ts)
# nobs is the number of observations
# m is the number of series
nobs, m = ts.shape
# Obtain the cointegrating vectors and corresponding eigenvalues
eigenvectors, eigenvalues = maximum_likelihood_estimation(ts, lags)
# Build table of critical values for the hypothesis test
critical_values_string = """2.71 3.84 6.63
13.43 15.50 19.94
27.07 29.80 35.46
44.49 47.86 54.68
65.82 69.82 77.82
91.11 95.75 104.96
120.37 125.61 135.97
153.63 159.53 171.09
190.88 197.37 210.06
232.11 239.25 253.24
277.38 285.14 300.29
326.53 334.98 351.25"""
select_critical_values = np.array(
critical_values_string.split(),
float).reshape(-1, 3)
critical_values = select_critical_values[:, 1]
# Calculate numbers of cointegrating relations for which
# the null hypothesis is rejected
rejected_r_values = []
for r in range(m):
if h_test(eigenvalues, r, nobs, lags, critical_values):
rejected_r_values.append(r)
return eigenvectors, rejected_r_values
def h_test(eigenvalues, r, nobs, lags, critical_values):
"""
Helper to execute the hypothesis test
"""
# Calculate statistic
t = nobs - lags - 1
m = len(eigenvalues)
statistic = -t * np.sum(np.log(np.ones(m) - eigenvalues)[r:])
# Get the critical value
critical_value = critical_values[m - r - 1]
# Finally, check the hypothesis
if statistic > critical_value:
return True
else:
return False
def maximum_likelihood_estimation(ts, lags):
"""
Obtain the cointegrating vectors and corresponding eigenvalues
"""
# Make sure we are working with array, convert if necessary
ts = np.asarray(ts)
# Calculate the differences of ts
ts_diff = np.diff(ts, axis=0)
# Calculate lags of ts_diff.
ts_diff_lags = lagmat(ts_diff, lags, trim='both')
# First lag of ts
ts_lag = lagmat(ts, 1, trim='both')
# Trim ts_diff and ts_lag
ts_diff = ts_diff[lags:]
ts_lag = ts_lag[lags:]
# Include intercept in the regressions
ones = np.ones((ts_diff_lags.shape[0], 1))
ts_diff_lags = np.append(ts_diff_lags, ones, axis=1)
# Calculate the residuals of the regressions of diffs and lags
# into ts_diff_lags
inverse = np.linalg.pinv(ts_diff_lags)
u = ts_diff - np.dot(ts_diff_lags, np.dot(inverse, ts_diff))
v = ts_lag - np.dot(ts_diff_lags, np.dot(inverse, ts_lag))
# Covariance matrices of the calculated residuals
t = ts_diff_lags.shape[0]
Svv = np.dot(v.T, v) / t
Suu = np.dot(u.T, u) / t
Suv = np.dot(u.T, v) / t
Svu = Suv.T
# ToDo: check for singular matrices and exit
Svv_inv = np.linalg.inv(Svv)
Suu_inv = np.linalg.inv(Suu)
# Calculate eigenvalues and eigenvectors of the product of covariances
cov_prod = np.dot(Svv_inv, np.dot(Svu, np.dot(Suu_inv, Suv)))
eigenvalues, eigenvectors = np.linalg.eig(cov_prod)
# Use Cholesky decomposition on eigenvectors
evec_Svv_evec = np.dot(eigenvectors.T, np.dot(Svv, eigenvectors))
cholesky_factor = np.linalg.cholesky(evec_Svv_evec)
eigenvectors = np.dot(eigenvectors, np.linalg.inv(cholesky_factor.T))
# Order the eigenvalues and eigenvectors
indices_ordered = np.argsort(eigenvalues)
indices_ordered = np.flipud(indices_ordered)
# Return the calculated values
return eigenvalues[indices_ordered], eigenvectors[:, indices_ordered]
johansen(df[['Brent','Gold','PMItot','PMIcomp','USD-EUR','SP','GDP']],3)